#include "StdAfx.h"
#include "SOM.h"
#include "time.h"
#include <complex>
#include "stdlib.h"
#include "math.h"
#include <iostream>
#include <fstream>
#include <atlimage.h>
#include "matrix.h"

#define natureNum 2.7182818284590452353602874713527
#define _CRT_RAND_S


//construct a SOM class
CSOM::CSOM(int row, int colm, int inputD, int maxC)
{
//initilized parameters
	srand((unsigned)time(NULL));
	sRow = row;
	sCol = colm;
	inputDim = inputD;
	outLayer = (som_node**)malloc(sizeof(som_node*)*row);
	for (int i=0; i<row; i++)
	{
		outLayer[i] = (som_node*)malloc(sizeof(som_node)*colm);
	}
	IniNode();

	is_loaded=0;
	currentIte = 0;
	maxCycle = maxC;
	mapRadius = max(row, colm)/4;
	train_destroy = 1;
	lprect = new RECT;
	is_illustrate = 0;
	inputSize = 0;
//determine the parameters
	iniNeighSize = mapRadius;
	iniLearningRate = 1;
	finalNeighSize = 1;
	finalLearningRate = 0.01;
	neighScale = log(iniNeighSize/finalNeighSize);
	learningScale = log(iniLearningRate/finalLearningRate);
//create an view of the SOM
	som_image = new CSOMImage(NULL);
	is_view = 1;
	if(!som_image->Create())
	{
		::AfxMessageBox(CString("Failed to create a view for the SOM"), MB_OK, 0);
		is_view = 0;
	}
	som_image->ShowWindow(SW_SHOW);
	CSOM::ViewCreate();
	cur_inputData = (double*)malloc(sizeof(double)*inputDim);
	Lattice = (double**)malloc(sizeof(double*)*sRow*sCol);
	int count = 0;
	for(int i=0; i<sCol; i++)
		for(int j=0; j<sRow; j++)
			Lattice[count++]=outLayer[j][i].weight;
	LatticePattern = new int[sRow*sCol];
}

//construct a SOM class
CSOM::CSOM(int row, int colm, int inputD, int maxC, double* maxInput, double* minInput)
{
//initilized parameters
	is_loaded=0;
	srand((unsigned)time(NULL));
	sRow = row;
	sCol = colm;
	inputDim = inputD;
	outLayer = (som_node**)malloc(sizeof(som_node*)*row);
	for (int i=0; i<row; i++)
	{
		outLayer[i] = (som_node*)malloc(sizeof(som_node)*colm);
	}
	IniNode(maxInput, minInput);

	currentIte = 0;
	maxCycle = maxC;
	mapRadius = max(row, colm)/4;
	train_destroy = 1;
	lprect = new RECT;
	is_illustrate = 0;
	inputSize = 0;
//determine the parameters
	iniNeighSize = mapRadius;
	iniLearningRate = 1;
	finalNeighSize = 1;
	finalLearningRate = 0.01;
	neighScale = log(iniNeighSize/finalNeighSize);
	learningScale = log(iniLearningRate/finalLearningRate);
//create an view of the SOM
	som_image = new CSOMImage(NULL);
	is_view = 1;
	if(!som_image->Create())
	{
		::AfxMessageBox(CString("Failed to create a view for the SOM"), MB_OK, 0);
		is_view = 0;
	}
	som_image->ShowWindow(SW_SHOW);
	CSOM::ViewCreate();
	cur_inputData = (double*)malloc(sizeof(double)*inputDim);
	Lattice = (double**)malloc(sizeof(double*)*sRow*sCol);
	int count = 0;
	for(int i=0; i<sCol; i++)
		for(int j=0; j<sRow; j++)
			Lattice[count++]=outLayer[j][i].weight;
	LatticePattern = new int[sRow*sCol];
}

//Load the SOM from a file
CSOM::CSOM(string filepath)
{
	is_loaded=1;
	//Load the data
	int counter = 0;
	try
	{
		//Create the map
		vector<vector<vector<double>>> somweights = bf.LoadBinaryData(filepath);
		sRow = (int)somweights.size();
		sCol = (int)somweights[0].size();
		inputDim = (int)somweights[0][0].size();
		outLayer = new som_node*[sRow];
		Lattice = new double*[sRow*sCol];
		for(unsigned int i=0; i<(unsigned int)sRow; i++)
		{
			outLayer[i] = new som_node[sCol];
			//Varify the columns
			if(somweights[i].size()!=sCol)
				throw string("Columns not correct");
			for(unsigned int j=0; j<(unsigned int)sCol; j++)
			{
				outLayer[i][j].X = i;
				outLayer[i][j].Y = j;
				unsigned int numberofdata;
				outLayer[i][j].weight = bf.VectorToPointer(somweights[i][j], numberofdata);
				Lattice[counter++] = outLayer[i][j].weight;
				//Varify the data dimension
				if(numberofdata!=inputDim)
					throw string("Data dimension not correct");
			}
		}
	}
	catch(string exceptionmessage)	//Catch the throw the exception message
	{
		throw exceptionmessage;
	}
}

//initial the nodes in SOM
void CSOM::IniNode(void)
{
	for (int i=0; i<sRow; i++)
	{
		for (int j=0; j<sCol; j++)
		{
			outLayer[i][j].X = i;
			outLayer[i][j].Y = j;
			outLayer[i][j].weight = ranDouble();
		}
	}
}

//initial the nodes in SOM
void CSOM::IniNode(double* maxvalue, double* minvalue)
{
	for (int i=0; i<sRow; i++)
	{
		for (int j=0; j<sCol; j++)
		{
			outLayer[i][j].X = i;
			outLayer[i][j].Y = j;
			outLayer[i][j].weight = ranDouble(maxvalue, minvalue);
		}
	}
}

//generate a random double array;
double* CSOM::ranDouble(void)
{
	double* ranNum = (double*)malloc(sizeof(double)*inputDim);
	for (int i=0; i<inputDim; i++)
	{
		try{
			ranNum[i] = double(abs(rand())+1)/double(abs(rand())+1);
			if(ranNum[i]>1)
				ranNum[i]=double(1)/ranNum[i];
		}catch(std::exception e){
			AfxMessageBox(CString(e.what()), MB_OK, 0);
		}
	}
	return ranNum;
}

//generate a random double array;
double* CSOM::ranDouble(double* maxvalue, double* minvalue)
{
	double* ranNum = (double*)malloc(sizeof(double)*inputDim);
	for (int i=0; i<inputDim; i++)
	{
		//Generate random doubles between 0~1
		ranNum[i] = double(abs(rand())+1.0)/double(abs(rand())+1.0);
		if(ranNum[i]>1.0)
			ranNum[i]=double(1)/ranNum[i];
		//Rescale the random doubles to minvalue~maxvalue
		ranNum[i] = ranNum[i]*maxvalue[i]+(1-ranNum[i])*minvalue[i];
	}
	return ranNum;
}

//computer distance between two double array.
double CSOM::wDistance(double* x, double* y)
{
	double dist = 0;
	for (int i=0; i<inputDim; i++)
	{
		dist += pow(x[i]-y[i], (double)2);
	}
	return dist;
}

//find the winner node
som_node* CSOM::findWinner(double* inputI)
{
	bool flag = 1;
	som_node* winner;
	double minDist = (double)ULONG_MAX;
	for (int i=0; i<sRow; i++)
	{
		for (int j=0; j<sCol; j++)
		{
			double temp = CSOM::wDistance(outLayer[i][j].weight, inputI);
			if (temp < minDist)
			{
				minDist = temp;
				winner = &outLayer[i][j];
				flag = 0;
			}
		}
	}
	//if the nodes are irregular or irregular input in the data, write to log and return error message
	if(flag)
	{
		::AfxMessageBox(CString("Can't find a winner for the input\r\nCheck the input data"), MB_OK, 0);
		int k;
		char* inputchar = bf.DoubleVectorToChar(inputI, inputDim, 10, ';', &k);
		bf.WriteToLog(CString(inputchar, k));
		CSOM::somWriteNodes(CString("som"));
		return &outLayer[0][0];
	}
	meanSquareError += minDist;
	return winner;
}

//Find the winners and store the distances
som_node* CSOM::findWinner(double* inputI, double& distance)
{
	bool flag = 1;
	som_node* winner;
	double minDist = (double)ULONG_MAX;
	for (int i=0; i<sRow; i++)
	{
		for (int j=0; j<sCol; j++)
		{
			double temp = CSOM::wDistance(outLayer[i][j].weight, inputI);
			if (temp < minDist)
			{
				minDist = temp;
				winner = &outLayer[i][j];
				flag = 0;
			}
		}
	}
	//if the nodes are irregular or irregular input in the data, write to log and return error message
	if(flag)
	{
		::AfxMessageBox(CString("Can't find a winner for the input\r\nCheck the input data"), MB_OK, 0);
		int k;
		char* inputchar = bf.DoubleVectorToChar(inputI, inputDim, 10, ';', &k);
		bf.WriteToLog(CString(inputchar, k));
		CSOM::somWriteNodes(CString("som"));
		return &outLayer[0][0];
	}
	meanSquareError += minDist;
	distance = minDist;
	return winner;
}

//update the weights
void CSOM::updateNode(som_node* node)
{
	int startx = (node->X-(int)neighsizedouble<0)?0:(node->X-(int)neighsizedouble);
	int endx = (node->X+(int)neighsizedouble>=sRow)?sRow:(node->X+(int)neighsizedouble+1);
	int starty = (node->Y-(int)neighsizedouble<0)?0:(node->Y-(int)neighsizedouble);
	int endy = (node->Y+(int)neighsizedouble>=sCol)?sCol:(node->Y+(int)neighsizedouble+1);
	//update the nodes and the colour
	if(is_illustrate==1)
	{
		for (int i=startx; i<endx; i++)
		{
			for (int j=starty; j<endy; j++)
			{
			//compare the distance
				double distTemp = pow((double)(node->X-outLayer[i][j].X), 2.0) + pow((double)(node->Y-outLayer[i][j].Y), 2.0);
				if (distTemp<double(neighsizedouble*neighsizedouble))
				{
					neighFunction = pow(natureNum, -distTemp/(2.0*neighsizedouble*neighsizedouble));
					//update the weights
					for (int k=0; k<inputDim; k++)
					{
						double temp = learningRate * neighFunction * (inputData[k][currentInput] - outLayer[i][j].weight[k]);
						outLayer[i][j].weight[k] += temp;
					}
					//show the image of the map
					//calculate the distance between the node and its neighbours
					double nodecolour = 0;
					for(int k=0; k<inputDim; k++)
						nodecolour += outLayer[i][j].weight[k];
					nodecolour /= inputDim;
					//calculate the colour
					nodecolour = nodecolour*(double)(255*255*32);
					nodecolour += (double)5;
					rgb rgbcolour;
					rgbcolour.rgbBlue = (byte)((int)(nodecolour)/(256*256));
					rgbcolour.rgbGreen = (byte)((int)(nodecolour)%(256*256)/256);
					rgbcolour.rgbRed = (byte)((int)(nodecolour)%(256*256)%256);
					CSOM::DrawNodes(node_pos[j*sRow+i].left, node_pos[j*sRow+i].top,
						node_pos[j*sRow+i].right, node_pos[j*sRow+i].bottom, rgbcolour);
				}
			}
		}
	}
	//update the nodes without changing the colour
	else
	{
		for (int i=startx; i<endx; i++)
		{
			for (int j=starty; j<endy; j++)
			{
			//compare the distance
				double distTemp = pow((double)(node->X-outLayer[i][j].X), 2.0) + pow((double)(node->Y-outLayer[i][j].Y), 2.0);
				if (distTemp<double(neighsizedouble*neighsizedouble))
				{
					neighFunction = pow(natureNum, -distTemp/(neighsizedouble*neighsizedouble*2.0));
					//update the weights
					for (int k=0; k<inputDim; k++)
					{
						double temp = learningRate * neighFunction * (inputData[k][currentInput] - outLayer[i][j].weight[k]);
						outLayer[i][j].weight[k] += temp;
					}
				}
			}
		}
	}
}

//update the node with the pattern colour
void CSOM::updateNode(som_node* node, rgb nodecolour)
{
	int startx = (node->X-(int)neighsizedouble<0)?0:(node->X-(int)neighsizedouble);
	int endx = (node->X+(int)neighsizedouble>=sRow)?sRow:(node->X+(int)neighsizedouble+1);
	int starty = (node->Y-(int)neighsizedouble<0)?0:(node->Y-(int)neighsizedouble);
	int endy = (node->Y+(int)neighsizedouble>=sCol)?sCol:(node->Y+(int)neighsizedouble+1);
	//update the nodes and the colour
	if(is_illustrate==1)
	{
		for (int i=startx; i<endx; i++)
		{
			for (int j=starty; j<endy; j++)
			{
				//compare the distance
				double distTemp = pow((double)(node->X-outLayer[i][j].X), 2.0) + pow((double)(node->Y-outLayer[i][j].Y), 2.0);
				if (distTemp<double(neighsizedouble*neighsizedouble))
				{
					neighFunction = pow(natureNum, -distTemp/(2.0*neighsizedouble*neighsizedouble));
					//update the weights
					for (int k=0; k<inputDim; k++)
					{
						double temp = learningRate * neighFunction * (inputData[k][currentInput] - outLayer[i][j].weight[k]);
						outLayer[i][j].weight[k] += temp;
					}
					//show the image of the map
					CSOM::DrawNodes(node_pos[j*sRow+i].left, node_pos[j*sRow+i].top,
						node_pos[j*sRow+i].right, node_pos[j*sRow+i].bottom, nodecolour);
				}
			}
		}
	}
	else
	{
		for (int i=startx; i<endx; i++)
		{
			for (int j=starty; j<endy; j++)
			{
				//compare the distance
				double distTemp = pow((double)(node->X-outLayer[i][j].X), 2.0) + pow((double)(node->Y-outLayer[i][j].Y), 2.0);
				if (distTemp<double(neighsizedouble*neighsizedouble))
				{
					neighFunction = pow(natureNum, -distTemp/(2.0*neighsizedouble*neighsizedouble));
					//update the weights
					for (int k=0; k<inputDim; k++)
					{
						double temp = learningRate * neighFunction * (inputData[k][currentInput] - outLayer[i][j].weight[k]);
						outLayer[i][j].weight[k] += temp;
					}
				}
			}
		}
	}
}

//training input data
void CSOM::somTraining(double** input, int inputS)
{
	meanSquareError = 0;
//if the SOM is loaded without training reject to train the SOM
	if(is_loaded==1)
	{
		::AfxMessageBox(CString("The SOM is loaded without training, can't be trained!!"), MB_OK, 0);
		return;
	}
//if previous train results not destroyed, destroy it
	if(inputS!=inputSize)
	{
		if(train_destroy==0)
			CSOM::trainDestroy();
		winners = (double**)malloc(sizeof(double*)*inputDim);
		for(int i=0; i<inputDim; i++)
		{
			winners[i] = (double*)malloc(sizeof(double)*inputS);
		}
		inputPat = (int*)malloc(sizeof(int)*inputS);
	}
//assign current input to som
	currentInput = 0;
	inputData = input;
	inputSize = inputS;
//construct the training results
//determine the number of iteration for calculate the neighing function and scaling
	numberIte = inputSize*maxCycle;

	trainingdistance.resize(inputS);
	while(currentInput<inputSize)
	{
		//change neighbouring function and learning scale
		learningRate = iniLearningRate * pow(natureNum, -(double)(currentIte)/(double)numberIte*learningScale);
		neighsizedouble = iniNeighSize * pow(natureNum, -(double)(currentIte)/(double)numberIte*neighScale);
		for(int i=0; i<inputDim; i++)
		{
			cur_inputData[i] = inputData[i][currentInput];
		}
		//find the winner
		winnerNode = CSOM::findWinner(cur_inputData, trainingdistance[currentInput]);
		//update weights
		CSOM::updateNode(winnerNode);
		//get the winners
		for(int i=0; i<inputDim; i++)
			winners[i][currentInput] = winnerNode->weight[i];
		//get the winner for the current input
		inputPat[currentInput] = winnerNode->X+winnerNode->Y*sRow;
		currentIte ++;
		currentInput ++;
		////Write the topology to file
		//if(currentIte%100==0)
		//	bf.AppendToFile(CString("C:\\somtopology.txt"), bf.douTostr(CSOM::somTopology())+CString("\t"));
	}
	//bf.AppendToFile(CString("C:\\somtopology.txt"), CString("\r\n"));
	train_destroy = 0;
	meanSquareError = meanSquareError / inputSize;
	return;
}

//get the winner at every iteration
double** CSOM::getTrainingWinners(void)
{
	if(train_destroy!=0)
	{
		::AfxMessageBox(CString("Erros!\r\nTraining results are destroyed"), MB_OK, 0);
		double** temp = (double**)malloc(sizeof(double*)*1);
		temp[0] = (double*)malloc(sizeof(double)*1);
		temp[0][0] = -1;
		return temp;
	}
	return winners;
}

//
CSOM::~CSOM(void)
{
}

//get the testing results
int CSOM::getTestingResults(double* testInput)
{
	int testResult;
	som_node* winnerNode = CSOM::findWinner(testInput);
	testResult = winnerNode->X + winnerNode->Y*sRow;
	return testResult;
}

int CSOM::getTestingResults(double *testInput, rgb testpattern)
{
	int testResult;
	som_node* winnerNode = CSOM::findWinner(testInput);
	testResult = winnerNode->X + winnerNode->Y*sRow;
	CSOM::DrawNodes(node_pos[testResult].left, node_pos[testResult].top, node_pos[testResult].right, node_pos[testResult].bottom, testpattern);
	return testResult;
}

//get the winner weight of the testing results
double* CSOM::getTestWinners(double* testInput)
{
	winnerNode = CSOM::findWinner(testInput);
	return winnerNode->weight;
}

//get the winner wieght of the testing results
double* CSOM::getTestWinners(double* testInput, rgb nodecolour)
{
	winnerNode = CSOM::findWinner(testInput);
	CSOM::DrawNodes(node_pos[winnerNode->X+winnerNode->Y*sRow].left, node_pos[winnerNode->X+winnerNode->Y*sRow].top,
		node_pos[winnerNode->X+winnerNode->Y*sRow].right, node_pos[winnerNode->X+winnerNode->Y*sRow].bottom, nodecolour);
	return winnerNode->weight;
}

//get the training results
int* CSOM::getTrainingResults()
{
	return inputPat;
}

//get the lattice
double** CSOM::getLattice(void)
{
	/*!!!!!!!!!!!!!!!!!! Remerber to transpose for the use of training !!!!!!!!!!!!!!!!!!!!!!!!!*/
    return Lattice;
}

//delete the som and clear the memory
void CSOM::somDestroy(void)
{
	if(train_destroy==0)
	{
		CSOM::trainDestroy();
	}
//destroy the output layer
	for(int i=0; i<sRow; i++)
	{
		for(int j=0; j<sCol; j++)
		{
			free(outLayer[i][j].weight);
		}
		free(outLayer[i]);
	}
	free(outLayer);
//if not illustrated, destroy the dialog
	if(is_illustrate == 0)
	{
		som_image->imageDestroy();
	}

	free(lpoha);
	free(cur_inputData);
	free(Lattice);
	free(LatticePattern);
	free(lprect);
	free(node_pos);
}

//destroy the training results
void CSOM::trainDestroy()
{
	if(train_destroy==1)
	{
		return;
	}
	for(int i=0; i<inputDim; i++)
	{
		free(winners[i]);
	}
	free(winners);
	free(inputPat);
	train_destroy = 1;
	return;
}

//return the input that's mapped to the node
int* CSOM::GetMapInput(int somsequence, int* mapnum)
{
//check if training results destroyed
	if(train_destroy==1)
	{
		int* temp = new int;
		*temp = -1;
		::AfxMessageBox(CString("ERROR!\r\nTraining Results Already Destroyed"), MB_OK, 0);
		return temp;
	}
//count the number of input mapped to the node
	int count = 0;
	for(int i=0; i<inputSize; i++)
	{
		if(inputPat[i] == somsequence)
		{
			count ++;
		}
	}
//construct the int pointer
	int* mapinput = (int*)malloc(sizeof(int)*count);
	count = 0;
	for(int i=0; i<inputSize; i++)
	{
		if(inputPat[i] == somsequence)
		{
			mapinput[count] = i;
			count ++;
		}
	}
	*mapnum = count;
	return mapinput;
}

//create the view for SOM
bool CSOM::ViewCreate()
{
//get the client size of the dialog;
	som_image->GetClientRect(lprect);
//create an image som for the dialog;
	if(!som_image->c_image.Create(lprect->right, lprect->bottom, 24, 0))
	{
		::AfxMessageBox(CString("failed to create view image"), MB_OK, 0);
		return 0;
	}
//initialized the image
	m_width = som_image->c_image.GetWidth();
	m_height = som_image->c_image.GetHeight();
	lpoha = (rgb**)malloc(sizeof(rgb*)*m_height);
	for(int j=0; j<m_height; j++)
	{
		lpoha[j] = (rgb*)som_image->c_image.GetPixelAddress(0, j);
	}
//plot the background
	for(int i=0; i<m_width; i++)
	{
		for(int j=0; j<m_height; j++)
		{
			lpoha[j][i].rgbRed = 103;
			lpoha[j][i].rgbGreen = 203;
			lpoha[j][i].rgbBlue = 211;
		}
	}
//calcualte the units of nodes
	u_x = m_width/(sCol+1);
	u_y = m_height/(sRow+1);
//plot the lines connecting the nodes
	for(int i=u_x; i<m_width-(m_width%(sCol+1)); i=i+u_x)
	{
		for(int j=u_y; j<m_height-(m_height%(sRow+1)+u_y); j++)
		{
			lpoha[j][i].rgbRed = 199;
			lpoha[j][i].rgbGreen = 230;
			lpoha[j][i].rgbBlue = 236;
		}
	}
	for(int j=u_y; j<m_height-(m_height%(sRow+1)); j=j+u_y)
	{
		for(int i=u_x; i<m_width-(m_width%(sCol+1)+u_x); i++)
		{
			lpoha[j][i].rgbRed = 199;
			lpoha[j][i].rgbGreen = 230;
			lpoha[j][i].rgbBlue = 236;
		}
	}
//get the position of nodes on the SOM
	rgb nodecolour;
	nodecolour.rgbRed = 113;
	nodecolour.rgbGreen = 213;
	nodecolour.rgbBlue = 211;
	node_pos = (RECT*)malloc(sizeof(RECT)*sRow*sCol);
	for(int i=0; i<sRow; i++)
	{
		for(int j=0; j<sCol; j++)
		{
			node_pos[sRow*j+i].left = u_x*(j+1) - u_x/4;
			node_pos[sRow*j+i].right = u_x*(j+1) + u_x/4;
			node_pos[sRow*j+i].top = u_y*(i+1) - u_y/4;
			node_pos[sRow*j+i].bottom = u_y*(i+1) + u_y/4;
		}
	}
//plot the nodes
	for(int i=0; i<sRow; i++)
	{
		for(int j=0; j<sCol; j++)
		{
			CSOM::DrawNodes(node_pos[j*sRow+i].left, node_pos[j*sRow+i].top,
				node_pos[j*sRow+i].right, node_pos[j*sRow+i].bottom, nodecolour);
		}
	}
	som_image->Invalidate();
	som_image->UpdateWindow();
	return 1;
}

//update the nodes colour
void CSOM::DrawNodes(int s_x, int s_y, int e_x, int e_y, rgb rgbcolour)
{
	for(int i=s_x; i<=e_x; i++)
	{
		for(int j=s_y; j<=e_y; j++)
		{
			lpoha[j][i] = rgbcolour;
		}
	}
	lprect->left = s_x;
	lprect->right = e_x+1;
	lprect->top = s_y;
	lprect->bottom = e_y+1;
	som_image->InvalidateRect(lprect, 1);
	som_image->UpdateWindow();
	return;
}

//move the image position
void CSOM::MoveImage(int pos)
{
	som_image->GetWindowRect(lprect);
	int w_w = lprect->right - lprect->left;
	int w_h = lprect->bottom - lprect->top;
	switch (pos)
	{
	case 1:		//move to the top left
		som_image->MoveWindow(0, 0, w_w, w_h, 1);
		break;
	case 2:		//move to the top middle
		som_image->MoveWindow(640-w_w/2, 0, w_w, w_h, 1);
		break;
	case 3:		//move to the top right
		som_image->MoveWindow(1280-w_w, 0, w_w, w_h, 1);
		break;
	case 4:		//move to the middle left
		som_image->MoveWindow(0, 512-w_h/2, w_w, w_h, 1);
		break;
	case 5:		//move to the middle right
		som_image->MoveWindow(1280-w_w, 512-w_h/2, w_w, w_h, 1);
		break;
	case 6:		//move to the bottom left
		som_image->MoveWindow(0, 1024-w_h, w_w, w_h, 1);
		break;
	case 7:		//move to the bottom middle
		som_image->MoveWindow(640-w_w/2, 1024-w_h, w_w, w_h, 1);
		break;
	case 8:		//move to the bottom right
		som_image->MoveWindow(1280-w_w, 1024-w_h, w_w, w_h, 1);
		break;
	default:	//move to the middle middle
		som_image->MoveWindow(640-w_w/2, 512-w_h/2, w_w, w_h, 1);
		break;
	}
	som_image->Invalidate();
	som_image->UpdateWindow();
}

//normalize the weight of node
void CSOM::normalizeWeight(void)
{
	double* weightAve = (double*)malloc(sizeof(double)*inputDim);
	double* weightSum = (double*)calloc(inputDim, sizeof(double));
	double* weightStd = (double*)malloc(sizeof(double)*inputDim);
	double nodeNum = sRow * sCol;
	for(int i=0; i<sRow; i++)
	{
		for(int j=0; j<sCol; j++)
		{
			for(int k=0; k<inputDim; k++)
			{
				weightSum[k] += outLayer[i][j].weight[k];
			}
		}
	}
	for(int k=0; k<inputDim; k++)
	{
		weightAve[k] = weightSum[k]/nodeNum;
		weightSum[k] = 0;
	}
	for(int i=0; i<sRow; i++)
	{
		for(int j=0; j<sCol; j++)
		{
			for(int k=0; k<inputDim; k++)
			{
				weightSum[k] += pow(outLayer[i][j].weight[k]-weightAve[k], 2);
			}
		}
	}
	for(int k=0; k<inputDim; k++)
	{
		weightStd[k] = sqrt(weightSum[k]/(nodeNum-1));
	}
	for(int i=0; i<sRow; i++)
	{
		for(int j=0; j<sCol; j++)
		{
			for(int k=0; k<inputDim; k++)
			{
				if(weightStd[k]!=0)
					outLayer[i][j].weight[k] = (outLayer[i][j].weight[k]-weightAve[k])/weightStd[k];
			}
		}
	}
	free(weightSum);
	free(weightAve);
	free(weightStd);
	return;
}

//write the output layer to file
void CSOM::somWriteNodes(CString str)
{
	//Header
	CString somstring = bf.douTostr(sRow*sCol);
	somstring += (CString(";")+bf.douTostr(inputDim));
	somstring += CString(";\r\n");
	//Body
	int sizeofchar;
	char* somchar = bf.DoubleMatrixToChar(Lattice, sRow*sCol, inputDim, 10, '\t', &sizeofchar);
	somstring += CString(somchar, sizeofchar);
	bf.stringToFile(somstring, str);
	free(somchar);
	return;
}

//Write the output layer to a file
void CSOM::somWriteToFile(CString filename)
{
	int sizeofchar;
	char* somstring = bf.DoubleMatrixToChar(Lattice, sRow*sCol, inputDim, 10, '\t', &sizeofchar);
	bf.WriteToFile(filename, CString(somstring, sizeofchar));
	return;
}

void CSOM::somSaveWeights(std::string filename)
{
	//Create the vectors
	vector<vector<vector<double>>> somweights;
	//Save the weights into the vectors
	somweights.resize(sRow);
	for(int i=0; i<sRow; i++)
	{
		somweights[i].resize(sCol);
		for(int j=0; j<sCol; j++)
		{
			somweights[i][j] = bf.PointerToVector(outLayer[i][j].weight, (unsigned int)inputDim);
		}
	}
	bf.SaveDataAsBinary(somweights, filename);
}

//determine if the map is illustrated
void CSOM::somIllustration(bool flag)
{
	is_illustrate = flag;
	return;
}

//load the SOM map from a double matrix
void CSOM::somLoadFromMatrix(double** m_Double, int m_SizeX, int m_SizeY, int m_DataDimension)
{
	//load the size and data dimension
	sRow = m_SizeX;
	sCol = m_SizeY;
	inputDim = m_DataDimension;
	CSOM::IniNode();
	//load the weights
	for(int i=0; i<sRow; i++)
	{
		for(int j=0; j<sCol; j++)
		{
			for(int k=0; k<m_DataDimension; k++)
			{
				outLayer[i][j].weight[k] = m_Double[i*sRow+j][k];
			}
		}
	}
	is_loaded = 1;
	return;
}

//Calcualte the topology of the map using Mu-Chun Su's method
double CSOM::somTopology(void)
{
	double mcsvalue = 0.0;
//	SOMCell[][] cells = som.getLnkSOMCell();
	for (int i = 0; i < sRow; i++)
	{
		for (int j = 0; j < sCol; j++)
		{
			int distance = 1;
			int left = j - distance;
			int right = j + distance;
			int top = i - distance;
			int bottom = i + distance;
			double* myweights = outLayer[i][j].weight;//cells[i][j].getWeights();
			double mypos[2] = {double(i), double(j)};
			//	go through all of the neurons that are outside of the immediate neighbors 
			// This is more than quartic complexity now, so there must be a better way to do it.  
			for (int rowout = 0; rowout < sRow; rowout++)
			{
				for (int colout = 0; colout < sCol; colout++)
				{
					// check to see if this element is not in the "outside" set
					if (rowout >= top && rowout <= bottom && colout >= left && colout <= right)
					{
						// so we are inside now, skip it!
						continue;
					}
					double outsidepos[2] = {double(rowout), double(colout)}; // stores the position of the outside cell
					// Now go through all the neighbors of the current cell and find the closest one to the outside cell
					double mingriddistanceout = double(ULONG_MAX);
					int winnerrow = 0;
					int winnercol = 0;
					for (int row = top; row <= bottom; row++)
					{
						if (row < 0 || row >= sRow)
								continue;
						for (int col = left; col <= right; col++)
						{
							if (col < 0 || col >= sCol)
								continue;
							if (row == i && col == j) // don't count yourself...
								continue;

							double neighborpos[2] = {double(row), double(col)};// stores the neighbors position
							// Find the distance between this neighbor and the outside cell
							double neighborgriddistance = bf.eDistance(&neighborpos[0], &outsidepos[0], 2);
							if (neighborgriddistance < mingriddistanceout)
							{
								mingriddistanceout = neighborgriddistance;
								winnerrow = row;
								winnercol = col;
							}
						}
					}

					// Now we have the closest neighbor cell, compute it's weights' distance to the outside cell's weights
					double* outsideweights = outLayer[rowout][colout].weight;
					double* neighweights = outLayer[winnerrow][winnercol].weight;

					double myweightdistout = CSOM::wDistance(myweights, outsideweights);
					double myweightdistneigh = CSOM::wDistance(myweights, neighweights);
					if (myweightdistout < myweightdistneigh)
					{ // add to the metric!!
						double mygriddistout = bf.eDistance(&mypos[0], &outsidepos[0], 2);
						//double scale = 1.0 - Math.exp(-Math.pow(mygriddistout, 2.0));
						double scale = 1.0 - pow(natureNum, (-pow(mygriddistout, 2.0)));
						double numerator = myweightdistneigh - myweightdistout;
						double denominator = myweightdistneigh;
				
						double value = scale * (numerator / denominator);
						mcsvalue += value;
					}
				}
			}
		}
	}
	return mcsvalue;
}

//Training the SOM with grouth truth data
void CSOM::somTraining(double **input, int inputS, int* groundtruth)
{
	meanSquareError = 0;
//if the SOM is loaded without training reject to train the SOM
	if(is_loaded==1)
	{
		::AfxMessageBox(CString("The SOM is loaded without training, can't be trained!!"), MB_OK, 0);
		return;
	}
//if previous train results not destroyed, destroy it
	if(inputS!=inputSize)
	{
		if(train_destroy==0)
			CSOM::trainDestroy();
		winners = (double**)malloc(sizeof(double*)*inputDim);
		for(int i=0; i<inputDim; i++)
		{
			winners[i] = (double*)malloc(sizeof(double)*inputS);
		}
		inputPat = (int*)malloc(sizeof(int)*inputS);
	}
//assign current input to som
	currentInput = 0;
	inputData = input;
	inputSize = inputS;
//construct the training results
//determine the number of iteration for calculate the neighing function and scaling
	numberIte = inputSize*maxCycle;
	//Get the number of patterns in the training
	int m_maxpattern = 1;
	for(int i=0; i<inputS; i++)
		if(m_maxpattern<groundtruth[i])
			m_maxpattern = groundtruth[i];
	if(m_maxpattern<=1)
	{
		::AfxMessageBox(CString("No patterns here"), MB_OK, 0);
		//return;
	}
	if(m_maxpattern>255)
	{
		::AfxMessageBox(CString("Too many patterns"), MB_OK, 0);
		return;
	}
	//Create the pattern colour
	rgb* patterncolour = new rgb[m_maxpattern];
	for(int i=0; i<m_maxpattern; i++)
	{
		patterncolour[i].rgbRed = (byte)((double)255/(double)(m_maxpattern-1)*(double)i);
		patterncolour[i].rgbGreen = (byte)((double)255/(double)(m_maxpattern-1)*(double)i);
		patterncolour[i].rgbBlue = (byte)((double)255/(double)(m_maxpattern-1)*(double)i);
	}
	//Create memory for the node pattern
	int** nodepattern = new int*[m_maxpattern];
	for(int i=0; i<m_maxpattern; i++)
		nodepattern[i] = (int*)calloc(sRow*sCol, sizeof(int));
//a pointer to current input data
	trainingdistance.resize(inputS);
	while(currentInput<inputSize)
	{
		//change neighbouring function and learning scale
		learningRate = iniLearningRate * pow(natureNum, -(double)(currentIte)/(double)numberIte*learningScale);
		neighsizedouble = iniNeighSize * pow(natureNum, -(double)(currentIte)/(double)numberIte*neighScale);
		for(int i=0; i<inputDim; i++)
		{
			cur_inputData[i] = inputData[i][currentInput];
		}
		//find the winner
		winnerNode = CSOM::findWinner(cur_inputData, trainingdistance[currentInput]);
		//update weights
		CSOM::updateNode(winnerNode, patterncolour[groundtruth[currentInput]-1]);
		//get the winners
		for(int i=0; i<inputDim; i++)
			winners[i][currentInput] = winnerNode->weight[i];
		//get the winner for the current input
		inputPat[currentInput] = winnerNode->X+winnerNode->Y*sRow;
		////Update the nodes
		//CSOM::DrawNodes(node_pos[inputPat[currentInput]].left, node_pos[inputPat[currentInput]].top,
		//	node_pos[inputPat[currentInput]].right, node_pos[inputPat[currentInput]].bottom,
		//	patterncolour[groundtruth[currentInput]-1]);
		//LatticePattern[inputPat[currentInput]]=groundtruth[currentInput];
		nodepattern[groundtruth[currentInput]-1][inputPat[currentInput]] ++;
		currentIte ++;
		currentInput ++;
		////Write the topology to file
		//if(currentIte%100==0)
		//	bf.AppendToFile(CString("C:\\somtopology.txt"), bf.douTostr(CSOM::somTopology())+CString("\t"));
	}
	//bf.AppendToFile(CString("C:\\somtopology.txt"), CString("\r\n"));
	//determine the pattern of each node and show the color
/*	for(int i=0; i<sRow*sCol; i++)
	{
		int currentpattern = -1;
		int maxpatternnumber = 0;
		for(int j=0; j<m_maxpattern; j++)
		{
			if(nodepattern[j][i]>maxpatternnumber)
			{
				maxpatternnumber = nodepattern[j][i];
				currentpattern = j;
			}
		}
		LatticePattern[i] = currentpattern;
	}
	for(int i=0; i<sRow*sCol; i++)
	{
		if(LatticePattern[i]>=0)
			continue;
		//if no winner has been mapped to the node, determine the node by its neighbors
		vector<int> temppatterns(m_maxpattern, -1);
//		for_each(temppatterns.begin(), temppatterns.end(), initializevalue<int>(-1));
		for(int k=i/sRow-1; k<=i/sRow+1; k++)
		{
			for(int j=i%sRow-1; j<=j%sRow+1; j++)
			{
				if(k>=0&&k<sRow&&j>=0&&j<sCol)
				{
					if(LatticePattern[k*sRow+j]>=0)
						temppatterns[LatticePattern[k*sRow+j]] ++;
				}
			}
		}
		LatticePattern[i] = static_cast<int>(distance(temppatterns.begin(), max_element(temppatterns.begin(), temppatterns.end())));
	}
	for(int i=0; i<sRow*sCol; i++)
	{
		//draw the nodes according to the color
		if(LatticePattern[i]>=0)
			CSOM::DrawNodes(node_pos[i].left, node_pos[i].top, node_pos[i].right, node_pos[i].bottom,
			patterncolour[LatticePattern[i]]);
		else
		{
			rgb rgbcolour = {255, 255, 255};
			CSOM::DrawNodes(node_pos[i].left, node_pos[i].top, node_pos[i].right, node_pos[i].bottom,
			rgbcolour);
		}
	}
	//Free memory and return;
*/	free(patterncolour);
	train_destroy = 0;
	meanSquareError = meanSquareError / inputSize;
	return;
}

//Return the patterns of the nodes on the SOM
int* CSOM::somGetLatticePattern(void)
{
	return LatticePattern;
}

//Save the image to files
void CSOM::somSaveImage(CString m_Destination)
{
	if(som_image!=NULL)
	{
		HRESULT hResult = som_image->c_image.Save(m_Destination);
		if(FAILED(hResult))
			bf.WriteToLog(CString("Can't save the image to file ")+m_Destination);
	}
	return;
}

//Get the mean square errors
double CSOM::somGetMeanSquareError(void)
{
	return meanSquareError;
}

//Get the multidimensional scaling metric
double CSOM::somGetMultidimensionalScaling(double** m_input, int m_inputsize)
{
	double* new_inputData = new double[inputDim];
	//Initialize the multidimensional scaling
	double multidimensional = 0;
	int cur_x;
	int cur_y;
	int new_x;
	int new_y;
	double distWinner;
	double distInput;
	//Loop to the find the winners and calculate earch multidimensional scaling
	for(int i=0; i<m_inputsize; i++)
	{
		//Save the current input vector to a double pointer
		for(int k=0; k<inputDim; k++)
			cur_inputData[k] = m_input[k][i];
		//Decompose the position of the current winner
		cur_x = inputPat[i]%sRow;
		cur_y = (inputPat[i]-cur_x)/sRow;
		for(int j=i+1; j<m_inputsize; j++)
		{
			//Save the new input vector to a double pointer
			for(int k=0; k<inputDim; k++)
				new_inputData[k] = m_input[k][j];
			//Decompose the position of the new winner
			new_x = inputPat[j]%sRow;
			new_y = (inputPat[j]-new_x)/sRow;
			//Compute the distance between the winners
			distWinner = ((double)cur_x-(double)new_x)*((double)cur_x-(double)new_x)+
				((double)cur_y-(double)new_y)*((double)cur_y-(double)new_y);
			//Compute the distance between the input
			distInput = CSOM::wDistance(new_inputData, cur_inputData);
			//Sum the multidimensional scaling
			multidimensional += pow(distInput-distWinner, (double)2);
		}
	}
	free(new_inputData);
	return multidimensional;
}

//Euclidean normalization
void CSOM::somEuclideanNormalization(void)
{
	for(int i=0; i<sRow; i++)
	{
		for(int j=0; j<sCol; j++)
		{
			//Calculate the square sum of the weights of each node
			double squaresum = 0;
			for(int k=0; k<inputDim; k++)
			{
				squaresum += outLayer[i][j].weight[k]*outLayer[i][j].weight[k];
			}
			if(squaresum == 0)
				squaresum = 1;
			else
				squaresum = sqrt(squaresum);
			//Normalize the weights of the node
			for(int k=0; k<inputDim; k++)
			{
				outLayer[i][j].weight[k] = outLayer[i][j].weight[k] / squaresum;
			}
		}
	}
	return;
}

//Partitioning the map
int* CSOM::somPartition(void)
{
/*----------------Calculate the threshold----------------*/
	int numberofpairs = 4*sRow*sCol-3*sRow-3*sCol+2;
	int kernelsizex = 3;
	int kernelsizey = 3;
	//Get the number of pairs
	double* distancevector = new double[numberofpairs];
	//Calculate the distances
	int m_count = 0;
	for(int i=0; i<sRow; i++)
	{
		for(int j=1; j<sCol; j++)
		{
			distancevector[m_count++] = CSOM::wDistance(outLayer[i][j].weight,
				outLayer[i][j-1].weight);
		}
	}
	for(int i=1; i<sRow; i++)
	{
		for(int j=0; j<sCol; j++)
		{
			distancevector[m_count++] = CSOM::wDistance(outLayer[i][j].weight,
				outLayer[i-1][j].weight);
		}
	}
	for(int i=1; i<sRow; i++)
	{
		for(int j=1; j<sCol; j++)
		{
			distancevector[m_count++] = CSOM::wDistance(outLayer[i][j].weight,
				outLayer[i-1][j-1].weight);
		}
	}
	for(int i=0; i<sRow-1; i++)
	{
		for(int j=1; j<sCol; j++)
		{
			distancevector[m_count++] = CSOM::wDistance(outLayer[i][j].weight,
				outLayer[i+1][j-1].weight);
		}
	}

	//Define a histogram
	int maxhistogram = (int)sqrt((double)numberofpairs);
	double* somhistogram = (double*)calloc(maxhistogram, sizeof(double));
	double maxofvector, minofvector;
	bf.VectorMaxMin(distancevector, numberofpairs, &maxofvector, &minofvector);
	bf.VectorRescale(distancevector, numberofpairs, 0, maxhistogram-1);

	//Get thei histogram
	for(int i=0; i<numberofpairs; i++)
		somhistogram[(int)distancevector[i]] += 1.0;

	//Normalize the histogram
	double sumofhistogram = 0;
	for(int i=0; i<maxhistogram; i++)
		sumofhistogram += somhistogram[i];
	for(int i=0; i<maxhistogram; i++)
		somhistogram[i] = somhistogram[i]/sumofhistogram;

	//Get the threshold and convert to the actual hisotgram
	double threshold = bf.GetOtsuThreshold(somhistogram, maxhistogram);
	threshold = threshold/maxhistogram*(maxofvector-minofvector) + minofvector;
	//set the threshold to be about 10% of the total
	threshold = ((double)sRow/(double)numberofpairs)*(maxofvector-minofvector)+minofvector;
/*-------------End of Calculate the threshold-------------*/

/*-------------Graphic bisection-------------*
	//Get the adjacency matrix
	int numberofneighbors = 4;
	double** adjacencymatrix = new double*[4];
	for(int i=0; i<numberofneighbors; i++)
	{
		adjacencymatrix = new double[sRow*sCol];
		int count = 0;
		for(int j=0; j<sCol; j++)
		{
			for(int i=0; i<sRow; i++)
			{
				//If the neighbor does not exist, set the matrix entry to 0, otherwise calcualte the
				//Euclidean distance between the two nodes
				if(j-1<0)
					adjacencymatrix[0][count] = -1.0;
				else
					adjacencymatrix[0][count] = CSOM::wDistance(outLayer[i][j].weight, outLayer[i][j-1].weight);
				if(i+1>=sRow)
					adjacencymatrix[1][count] = -1.0;
				else
					adjacencymatrix[1][count] = CSOM::wDistance(outLayer[i][j].weight, outLayer[i+1][j].weight);
				if(j+1>=sCol)
					adjacencymatrix[2][count] = -1.0;
				else
					adjacencymatrix[2][count] = CSOM::wDistance(outLayer[i][j].weight, outLayer[i][j+1].weight);
				if(i-1<0)
					adjacencymatrix[3][count] = -1.0;
				else
					adjacencymatrix[3][count] = CSOM::wDistance(outLayer[i][j].weight, outLayer[i-1][j].weight);
				count ++;
			}
		}
	}

	//Define three measurement: association between A and V, association between B and V, association between A and B
	double oldassociationAV = 0.0;
	double oldassociationBV = 0.0;
	double oldassociationAB = 0.0;

	//Create status for the nodes A=0, B=1;
	int** nodeclass = new double*[sRow];
	for(int i=0; i<sRow; i++)
	{
		nodeclass[i] = (int*)calloc(sCol, sizeof(int));
	}

	//Initialize the associations: A=V, B=0
	for(int i=0; i<numberofneighbors; i++)
	{
		for(int j=0; j<sRow*sCol; i++)
		{
			if(adjacencymatrix[i][j]>=0)
			{
				associationAV += adjacencymatrix[i][j];
			}
		}
	}

	//Create association for a new state
	double newassociationAV = oldassociationAV;
	double newassociationBV = oldassociationBV;
	double newassociationAB = oldassociationAB;

	double oldcut = 1.0;
	double newcut = 1.0;
/*-------------Graphic bisection-------------*/

	//Create the pointers
	int* SOMPartition = (int*)calloc(sRow*sCol, sizeof(int));
	//Set the status of each node
	bool** NodeStatus = new bool*[sRow];
	for(int i=0; i<sRow; i++)
		NodeStatus[i] = (bool*)calloc(sCol, sizeof(bool));

	//Find the starting node with the smallest distance to its neighbors
	CPoint cpoint(0, 0);
	double mindistance = (double)ULONG_MAX;
	int count = 0;
	double tempdistance;
	for(int i=1; i<sRow-1; i++)
	{
		for(int j=1; j<sCol-1; j++)
		{
			//Get the sum distance between its neighbors
			count = 0;
			tempdistance = 0;
			for(int k=0; k<3; k++)
			{
				for(int l=0; l<3; l++)
				{
					tempdistance += CSOM::wDistance(outLayer[i][j].weight,
						outLayer[i+k-1][j+l-1].weight);
				}
			}
			//Determine the min distance
			tempdistance = tempdistance/9.0;
			if(mindistance>tempdistance)
			{
				mindistance = tempdistance;
				cpoint.x = i;
				cpoint.y = j;
			}
		}
	}

	//Create to vectors for merging
	vector<CPoint>Region1;
	vector<CPoint>Region2;
	Region1.push_back(cpoint);
	SOMPartition[cpoint.y*sRow+cpoint.x] = 1;
	NodeStatus[cpoint.x][cpoint.y] = 1;
	//nodeclass[cpoint.x][cpoint.y] = 1;
	//for(int i=0; i<numberofneighbors; i++)
	//{
	//	newassociationAV -= adjacencymatrix[i][cpoint.x+cpoint.y*sRow];
	//	newassociationAB += adjacencymatrix[i][cpoint.x+cpoint.y*sRow];
	//}

	//Create indicator
	bool flag = 1;
	while(flag)
	{
		//Pop data in Region1 to merge
		while(Region1.size()!=0)
		{
			CPoint CurrentPoint = Region1.back();
			Region1.pop_back();			
			////Check the status of merging the one on top
			//if(CurrentPoint.y-1>=0)
			//{
			//	oldassociationAV = newassociationAV;
			//	oldassociationBV = newassociationBV;
			//	oldassociationAB = newassociationAV;
			//	for(int i=0; i<numberofneighbors; i++)
			//	{
			//	}
			//}
			for(int i=0; i<3; i++)
			{
				for(int j=0; j<3; j++)
				{
					if((CurrentPoint.x+i-1)>=0&&(CurrentPoint.x+i-1)<sRow&&
						(CurrentPoint.y+j-1)>=0&&(CurrentPoint.y+j-1)<sCol)
					{
						//If a node to be merged, push it to Region2 to be used for merging
						if(NodeStatus[CurrentPoint.x+i-1][CurrentPoint.y+j-1]==0)
						{
							if(CSOM::wDistance(outLayer[CurrentPoint.x+i-1][CurrentPoint.y+j-1].weight,
								outLayer[CurrentPoint.x][CurrentPoint.y].weight)<threshold)
							{
								NodeStatus[CurrentPoint.x+i-1][CurrentPoint.y+j-1] = 1;
								SOMPartition[(CurrentPoint.x+i-1)+(CurrentPoint.y+j-1)*sRow] = 1;
								CPoint* NewPoint = new CPoint(CurrentPoint.x+i-1, CurrentPoint.y+j-1);
								Region2.push_back(*NewPoint);
							}
							else	//If not, set the node to be checked
							{
								NodeStatus[CurrentPoint.x+i-1][CurrentPoint.y+j-1] = 1;
								SOMPartition[(CurrentPoint.x+i-1)+(CurrentPoint.y+j-1)*sRow] = 2;
							}
						}
					}
				}
			}
		}
		//Pop data in Region1 to merge
		while(Region2.size()!=0)
		{
			CPoint CurrentPoint = Region2.back();
			Region2.pop_back();
			for(int i=0; i<3; i++)
			{
				for(int j=0; j<3; j++)
				{
					if((CurrentPoint.x+i-1)>=0&&(CurrentPoint.x+i-1)<sRow&&
						(CurrentPoint.y+j-1)>=0&&(CurrentPoint.y+j-1)<sCol)
					{
						//If a node to be merged, push it to Region1 to be used for merging
						if(NodeStatus[CurrentPoint.x+i-1][CurrentPoint.y+j-1]==0)
						{
							if(CSOM::wDistance(outLayer[CurrentPoint.x+i-1][CurrentPoint.y+j-1].weight,
								outLayer[CurrentPoint.x][CurrentPoint.y].weight)<threshold)
							{
								NodeStatus[CurrentPoint.x+i-1][CurrentPoint.y+j-1] = 1;
								SOMPartition[(CurrentPoint.x+i-1)+(CurrentPoint.y+j-1)*sRow] = 1;
								CPoint* NewPoint = new CPoint(CurrentPoint.x+i-1, CurrentPoint.y+j-1);
								Region1.push_back(*NewPoint);
							}
							else	//If not, set the node to be checked
							{
								NodeStatus[CurrentPoint.x+i-1][CurrentPoint.y+j-1] = 1;
								SOMPartition[(CurrentPoint.x+i-1)+(CurrentPoint.y+j-1)*sRow] = 2;
							}
						}
					}
				}
			}
		}
		//Check the flag
		if(Region1.size()==0&&Region2.size()==0)
			flag = 0;
	}

	for(int i=0; i<sRow*sCol; i++)
	{
		if(SOMPartition[i]!=1)
			SOMPartition[i] = 2;
	}
	return SOMPartition;
}

//Show the partition of the map
void CSOM::somShowPartition(void)
{
	int* sompartition = CSOM::somPartition();
	for(int i=0; i<sRow*sCol; i++)
	{
		if(sompartition[i]==1)
		{
			rgb rgbcolour ={(byte)255, (byte)255, (byte)255};
			CSOM::DrawNodes(node_pos[i].left, node_pos[i].top, node_pos[i].right,
				node_pos[i].bottom, rgbcolour);
		}
		else
		{
			rgb rgbcolour = {(byte)0, (byte)0, (byte)0};
			CSOM::DrawNodes(node_pos[i].left, node_pos[i].top, node_pos[i].right,
				node_pos[i].bottom, rgbcolour);
		}
	}
	free(sompartition);
}

//Write the som into an image if the input dim is 3
void CSOM::somWriteToImage(CString filename)
{
	CImage somimage;
	somimage.Create(sCol, sRow, 24, 0);
	rgb** lpoha = new rgb*[sRow];
	for(int j=0; j<sCol; j++)
	{
		lpoha[j] = (rgb*)somimage.GetPixelAddress(0, j);
	}
	double maxnode, minnode;
	bf.MatrixMaxMin(Lattice, sRow*sCol, inputDim, &maxnode, &minnode);
	bf.MatrixRescale(Lattice, sRow*sCol, inputDim, 0.0, 255.0);
	for(int j=0; j<sRow; j++)
	{
		for(int i=0; i<sCol; i++)
		{
			lpoha[j][i].rgbRed = (byte)outLayer[j][i].weight[0];
			lpoha[j][i].rgbGreen = (byte)outLayer[j][i].weight[1];
			lpoha[j][i].rgbBlue = (byte)outLayer[j][i].weight[2];
		}
	}
	bf.MatrixRescale(Lattice, sRow*sCol, inputDim, minnode, maxnode);
	somimage.Save(filename);
}

//Plot the index of the SOM
void CSOM::somPlotIndex(int* somindex, int indexlength)
{
	//Verity the index length
	if(indexlength!=sRow*sCol)
	{
		::AfxMessageBox(CString("Index length not correct!"), MB_OK, 0);
		return;
	}
	
	//Let us plot the index with different grayscales
	int count = 0;
	for(int i=0; i<sCol; i++)
	{
		for(int j=0; j<sRow; j++)
		{
			rgb rgbcolour;
			rgbcolour.rgbRed = somindex[count];
			rgbcolour.rgbGreen = somindex[count];
			rgbcolour.rgbBlue = somindex[count];
			CSOM::DrawNodes(node_pos[count].left, node_pos[count].top, node_pos[count].right,
				node_pos[count].bottom, rgbcolour);
			count ++;
		}
	}
	
	return;
}

//Get the size of the SOM
void CSOM::somGetSize(unsigned int &rows, unsigned int &columns)
{
	rows = (unsigned int)CSOM::sRow;
	columns = (unsigned int)CSOM::sCol;
	return;
}

//Get training distance
vector<double> CSOM::somGetTrainingDistance()
{
	return trainingdistance;
}

//Get the size of the map
_2dpoint CSOM::GetMapSize()
{
	_2dpoint sizeofmap = {(unsigned int)sRow, (unsigned int)sCol};
	return sizeofmap;
}

/*force to change the learning rate and neighborhood function and training cycle of the SOM*/
/*careful to use this, it might cause the SOM to behave stupidly*/
void CSOM::ForceChangeLR(double learningrate)
{
	CSOM::currentIte = 0;
	CSOM::learningRate = learningrate;
	learningScale = log(iniLearningRate/finalLearningRate);
	return;
}

void CSOM::ForceChangeNF(int neighborhoodsize)
{
	CSOM::currentIte = 0;
	CSOM::neighSize = neighborhoodsize;
	neighScale = log(iniNeighSize/finalNeighSize);
	return;
}

void CSOM::ForceChangeTC(int trainingcycle)
{
	CSOM::currentIte = 0;
	CSOM::maxCycle = trainingcycle;
	return;
}
/*-----------------------------------------------------------------------------------------------*/